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Abstract 

In this paper we report on a Monte Carlo study of a diluted Ising antiferro- 
magnet on a fee lattice. This is a typical model example of a highly frustrated 
antiferromagnet, and we ask, whether sufficient random dilution of spins does 
produce a spin glass phase. Our data strongly indicate the existence of a spin 
glass transition for spin-concentration p < 0.75: We find a divergent spin 
glass susceptibility and a divergent spin glass correlation length, whereas the 
antiferromagnetic correlation length saturates in this regime. Furthermore, 
we find a first order phase transition to an antiferromagnet for 1 > p > 0.85, 
which becomes continuous in the range 0.85 > p > 0.75. Finite size scal- 
ing is employed to obtain critical exponents. We compare our results with 
experimental systems as diluted frustrated antiferromagnets as Zni_pMnpTe. 
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I. INTRODUCTION 



The family of diluted magnetic semiconductors (DMS) of the general form A^^_pMnpB^^ 
encompasses a wide variety of alloys which have been under extensive investigation during 
the past 15 years. These alloys form a zincblende structure, where the B^^ element occupies 
one fee lattice while A^^ and Mn share the second fee lattice. One fundamental aspect of re- 
search has focused on the magnetic order of these systems, since they offer practical examples 
of strongly frustrated, randomly diluted three dimensional fee Heisenberg antiferromagnets 
(afm) with dominant nearest neighbour interaction^'!. 

In this paper we present results of a Monte Carlo study of a diluted frustrated Ising 
model on a fee lattice given by the Hamiltonian 

{1 with prob. p 
(1) 
with prob. 1 — p. 

Here, J is the coupling constant, which we will set J = — 1 henceforward, and p G [0, 1] is 
the probability that a lattice site i is occupied with an Ising-spin Sj. We are interested in 
the static properties of this model for different dilution regimes. Besides the pure {p = 1) 
and the slightly diluted case (p ~ 1), that has already been studied by MC simulationiFi 
and other methodsi0, we concentrate our interest on the strong dilution regime, which has 
only been investigated in experimental Heisenberg-systems as mentioned above. Although 
in our work we perform a simulation of an Ising system we find that some typical DMS 
results can be reproduced with our simplified model. 

Neutron diffraction experiments of thin Zni_pMnpTe films for p G [1.0,0.85]) revealed 
a first order phase transition to an antiferromagnetically ordered state of "Type-Ill"!. At 
approximately p = 0.85 a tricritical point is encountered, where the phase transition becomes 
continuous. The antiferromagnetic order in the regime p G [0.75,0.85] is truly long range, 
however, below p 0.75 a transition to a short range ordered phase seems to appeaii'i. 

Below p = 0.7 most experimental results have led to the view, that one encounters a tran- 
sition to a spin-glass-like phase at a fairly well defined temperature T^. Characteristic spin 



glass features are: (i) remanence effects in the frozen state0'0, (ii) a pronounced cusp in tlie 
susceptibility around with strong frequency dependence of the cusp temperatureS, (iii) 
absence of long range spin order as observed by magnetic neutron diffractionl^, (iv) dynamic 
scaling near Tc of frequency-dependent response function^Uli and — most importantly — 
(v) a divergent nonlinear susceptibility around the cusp temperature^. 

On the other hand, the antiferromagnetic correlation length ^afm grows continuously 
with decreasing temperature until it saturates at the cusp temperature to an enormously 
large value as high as 70 A at p = 0.7; it is only below p = 0.4 that short range afm order 
disappears. In the intermediate dilution range p G (0.4, 0.7) the spin glass interpretation 
has been questioned, and it was suggested briefiy0 (on the basis of an "activated" scaling 
analysis) that the equilibrium transition at was to the antiferromagnetically ordered state 
of a random-field system. 

This motivates a numerical study of a diluted antiferromagnet, in which we can observe 
the interplay of strong afm local order with spin glass order, and can measure the quantities 
now considered to be the signatures for spin-glass transitions. Because Heisenberg model 
simulations demand more computer time, and because of the more convoluted controversies 
regarding the existence of a sharp phase transition for continuous spins, we adopt here an 
Ising model; in addition, we have retained only the first-neighbour exchange interaction. 

The paper is organized as follows: Section || is concerned with the main theoretical 
arguments that guide our expectations for the results of our simulations in the distinct 
dilution regimes. In particular, we discuss the possible universality class of the proposed 
spin glass phase. Section |T| describes technical aspects of our simulations, in particular 



the equilibration criterion. In section IV -VI of our report we shall present representative 



data for the distinct dilution regimes. In section ^ we concentrate on the pure case and 
on weak dilution (p ~ 1), where we investigate how the order of the transition is being 
modified by disorder in form of stochastically removing spins from the lattice. Section |V| is 
concerned with the regime p = 0.8, where a continuous phase transition with afm ordering 



is found; critical exponents are determined by finite size scaling. Section ^ investigates the 



intermediate and strong dilution range {p G [0.3,0.7]), where the question of the magnetic 
ordering is our main concern. A summary of the results, a comparison to experimental 
systems and our final conclusions will be given in section |V11| . 



II. THEORETICAL BACKGROUND 

In this section, we collect the qualitative expectations for all the different concentration 
regimes expected in the phase diagram. As with the results in the subsequent sections, we 
begin by reviewing the pure case and proceed in the direction of greater dilution. The global 



phase diagram is qualitatively similar to those conjectured for vector spins in Ref. |T8| and 
Ref. |19|, except of course that distinctive collinear and noncollinear phases cannot exist in 
the Ising case. 

A. Pure fee Ising antiferromagnet 

The pure Ising antiferromagnet on a fee lattice has been extensively studied both ana- 
lytically and with simulations. Each spin has 12 nearest neighbours which in the ground 
state can only satisfy 8 bonds, 4 of them being always violated. This effect of frustration, 
which follows from the triangles in the fee lattice, leads to a large ground state degeneracy^ 
of the order 0{2^), where L denotes the linear extent of the system. Thus the entropy per 
spin is zero as L — oo. 

At small temperatures in this system, thermal fluctuations generate free energy terms 
which have the same effect as ferromagnetic second neighbor interactions: this favors the 
"Type-I" afm order, meaning that the system orders into one of the six periodic ground 
states with a (100) type ordering wavevectori. This is an example of what Villain called 
"ordering due to disorder"iil. 

The discrete choice between the three (100) type directions suggests a similarity in be- 
havior to the 3-state Potts model, which has a weakly first order phase transition in three 
dimensions. Indeed, the 4 — e renormalization group predicts a first-order phase transitionil. 



Simulation^li'ii and series expansionsi confirm that the phase transition is at finite tem- 
perature and is of first order. 

The antiferromagnetic state may be handled quantitatively by constructing (in the spirit 
of Ref. |23) a three-component staggered-magnetization order parameter = {ml, ml, ml) 



with components 

K = ^ E exp(-i Tj ■ k^) Sj (2) 

j 

for /i = 1, 2, 3; here the ordering wave vectors are fci = (vr, 0, 0) and cyclic permutations (We 
have taken the lattice constant to be unity.). 



B. Weak dilution: random-field effects 

Dilution in frustrated systems (without any external field) couples to the order parameter 
as a random field does in a ferromagnetlllii. Take the case of rather weak dilution, which 
justifies assuming (as a sort of variational state) one of the six (100) type ground states. 
Consider the effect of strengthening one bond lying within the xy plane: it will favor the 
four states with (100) and (010) ordering wavevectors and disfavor the two with (001) 
wavevectors, since the bond in question is violated in the (001) states. The effect is much 
like a random field, except that it does not destroy the global up/down symmetry. In our 
case of site dilution, the random-field-like effects of removing an isolated site cancel each 
other; however, removing a pair of sites has the same effect as would strengthening the bond 
between themSi. 

Quite generally, when the random field is sufficiently strong, the first-order transition 
is converted to a continuous one0. In the present context, since the effective random field 
grows with dilution, this argument predicts a tricritical point: the ordering transition is 
first-order for p > ptn but becomes continuous for p < ptr^- 

For p < ptri the transition from the paramagnet is expected to be a novel universality 
clasl^. It would seem plausible if its dynamic scaling behavior were of the "activated" type, 



5 



as in the random-field Ising model. No experiments have tested this, however (The materials 
in this concentration range, roughly 0.7 < p < 0.85, can be grown only as thin epitaxial 
slabs, meaning that very little signal is available for susceptibility experiments.). A claim 
was made that "activated" scaling could fit the data0 for lower values of p, which we would 
identify as the spin-glass phase, but this was quickly corrected0@lii. 

When the "effective random fields" are sufficiently strong, the antiferromagnetic order 
disappears and is replaced by spin-glass order at a critical value p^. Note that this thresh- 
old to propagate afm order is far above the pc, the geometrical percolation threshold for 
propagating connectivity of nearest-neighbor sitesS (On the fee lattice, pc = 0. 1950.). 



C. Spin glass phase 

Any spin glass, by definition, requires random frustration. This can be realized by 
dilution of a uniformly frustrated antiferromagnet, as in the present case, just as well as by 
a random mix of ferromagnetic and antiferromagnetic interactions^. Indeed, the effective 
coupling between two spins may be ferromagnetic or antiferromagnetic depending on how 
intervening sites happen to be occupied. Of course, this spin glass state is expected to show 
residual short-range correlations as opposed to the cases of the ±J or Gaussian random 
bond distributions, where symmetry implies that [(sjSj)] = if z 7^ j. 

It is intriguing that this "spin glass" state, which occurs below p^,, from the viewpoint 
of the antiferromagnetic phase, is the same as the disordered domain state which is favored 
by the "effective random fields" mentioned aboveil. This state is different from the familiar 
random-field disordered state (and similar to the usual spin glass state) because the global 
up/down symmetry is preserved; consequently, ioi p < there is still a true phase transition 
in which this symmetry is brokenil. 
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D. Universality 



In a concentration range (p'^p*), the ground state is presumed to be a spin glass. Spin 
glass investigations tend still to be preoccupied with the issue of the existence of a transition 
as a function of dimension, external field, and spin type. Indeed, it is still unsettled whether 
the d = 3 Ising spin glass really has a transition at finite temperature, or whether it is at 
the lower critical dimensionalityS'ii. Rather little has been done to test the universality of 
the critical exponents, as almost all simulations have used simple cubic (sc) lattices with the 
discrete ± J distribution of random bonds. Monte Carlo and series studies for the fee lattice 
with ±J bondsH gave values of the spin-glass exponents r], z/, and 7 consistent with the 
sc ±J model; so did diluted ±J model on a simple cubic lattice!! and (modulo large error 
bars) Gaussian-distributed random bonds on the simple cubic latticed. The above results 
are consistent with universality. However, it has also been proposed that rj is more negative 
and the Binder cumulant is larger for Gaussian bond randomness than for ± J randomness^; 
presumably the diluted fee is more similar to the latter model, since its discrete randomness 
generically allows exact degeneracy of ground states. 



E. Theory of p'^ (spin glass near percolation) 

We now consider the approach p p'^, where the spin glass long-range order finally 
disappears. In this regime, the order is just barely propagating along tortuous, effectively 
one-dimensional paths, and consequently we expect Tc (exponentially) as p ^ p'^HS. 

Note that p'^ > p^. In frustrated models with discrete bond distributions, such as the 
present case, two portions of a connected cluster might be connected by (say) two chains of 
bonds, each canceling the other and allowing one portion to be flipped relative to the other 
portion at no cost in energy; for propagation of order, it is as if no chain existed, i. e., the 
effective concentration of bonds is lowered by the cancellations. 
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III. TECHNICALITIES 



We use the single spin flip Monte Carlo Metropolis algorithm in our simulations. Spins 
are updated sequentially and randomly. Periodic boundary conditions are imposed, limiting 
the possible lattice sizes to even numbers. Spins are represented on a cubic lattice with 
next nearest neighbour interactions to obtain a fee lattice. Therefore, every lattice of linear 
size L contains N = sites. We simulated lattice sizes L = 4,6,8,10 with M ^ 120 

realizations of the disorder and L = 16 with M = 40. We investigate the model in the 
concentration range p G [0.3, 1.0]. 

A standard criterion by Bhatt and Young0 was applied to test the equilibration of the 
systems throughout the whole simulation, where we observe a continuous phase transition. 
The procedure is to obtain two estimates of the spin glass phase indicator, the spin glass 
susceptibility 



Here and later, the brackets (■ ■ ■) denote thermal averaging and [■ ■ ■] configurational aver- 
aging. We obtain xsg by calculating the second moment of the spin glass order parameter 
defined in the two alternative ways, (i) as the overlap 



Here, s- and s\ denote two sets of spins (replicas) with the same set of occupied sites and 
uncorrelated initialization and to is the time initially used for equilibration. 
With these definitions, we can compute two estimates of xsg as follows, i. e.. 




(3) 




(4) 



and (ii) as the autocorrelation (self-overlap) 




(5) 





respectively the four-spin correlation function. 
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(b) 
XsG 



N 



(7) 



where (■ ■ ■)r = 7 EU" " and r = to- 

The equihbration time to was raised on a logarithmic scale and we only accepted a run, 
if both estimates of xsg agreed after this time within certain limits, typically of the order 
of 5% of their joint mean value. The longest runs performed were up to 2 * 10^ Monte Carlo 
steps per spin (MCS). Most of the simulations were performed on HP workstations at the 
Institut fiir Theoretische Physik, Gottingen and on the Intel-Paragon parallel computer at 
the Hochstleistungsrechenzentrum Jiilich. The program was parallelised by using PVM 3.2 
software^ in order to simulate many systems simultaneously. 



IV. FIRST ORDER ANTIFERROMAGNETIC TRANSITION (WEAK 

DILUTION) 

In this section we investigate the pure afm and the slightly diluted regime, i. e., p G 
[0.85,1.0]. We wish to determine the order of the phase transition and how the order of 
the transition is changed by introducing disorder into the system in form of slight stochastic 
dilution. The magnetic order in this regime is clearly antiferro magnet, consistent with earlier 
simulations, and shall be closer examined in section 0. 

The pure antiferromagnet on a fee lattice is known to undergo a temperature driven 
first order phase transitioni&il'ii, as mentioned in the previous section. Early Monte Carlo 
simulations by Crest and Cabi as well as Ciebultowic^ reported a change from a first 
order to a continuous phase transition upon dilution. Crest and Cabl located the tricritical 
point at a critical concentration ptri = 0.93 using Ising-spins, whereas Ciebultowicz found 
a slightly lower ptri = 0.85 with a Heisenberg-spin simulation. However, in both of those 
simulations, no averaging over the disorder was performed, so that we re-investigated this 
regime. 

An important quantity for a first order phase transition is the internal energy density 
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At the critical temperature, this quantity indicates a first order phase transition by a discon- 
tinuity (latent heat), which can be seen in Fig. |^, where [(e)] is plotted versus temperature 
T; with increasing lattice size a pronounced discontinuity at the transition temperature can 
be observed, revealing clearly a first order transition. 

In Fig. ^ we present our data of the internal energy density for p = 0.9. For this 
concentration the transition appears to be continuous. Apparently, a tricritical point is 
hard to locate with these Monte Carlo data, since it is not clear, whether in the limit of 
smaller temperature steps and larger systems the transition will turn out to be continuous 
or not. In order to determine the tricritical concentration more efficiently we used a method 
first introduced by Lee & Kosterlitz0, which is based on a histogram algorithm by Ferrenberg 
& SwendsenS. 



A. Histogram algorithm and finite-size scaling analysis method 

This method^ uses one long Monte Carlo run to estimate the free energy at several 
temperatures close to Tc. During the runs a histogram of the internal energy density e is 
accumulated . This histogram serves as an estimate of the equilibrium energy distribution 
(after correct normalization) 

P^(e) = ^exp(S(e)-/5e). (9) 

Here, Zf^ is the partition function at the inverse temperature f3 = l/ksT and S{e) is the 
entropy. Notice, that Ppi^e) (||) is proportional to exp(— /5F(e)), where F denotes the free en- 
ergy. The distribution Pg(e) can now be used to generate the distribution (and consequently 
F{e)) at a different inverse temperature (3' in the vicinity of (3. 

We are interested in the situation in which F{e) has two minima ei and 62 (i.e. (^ 
has maxima at these energies). The temperature at which F{ei,L) = F{e2,L) is taken 
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to be the effective critical temperature for the given size L. Then we evaluate the "gap" 
AF = F{em) — F{ei^2) between the free energy values at those two energies and at the 
maximum Cm in between them. From a scaling analysis of AF, one can identify whether a 
phase transition is first order@. 

A state with e G [ei, 62] consists of a domain of ordered and a domain of disordered phase 
coexisting, separated by a (c/ — l)-dimensional interface surrounding the droplet of minority 
phase. Therefore, one can expand the free energy 

Fie,L) = L'fo{e) + L''-'Me) + 0{L''-'), (10) 

where the bulk free energy density /o is minimal and constant for e G [ei, 62] and the surface 
term /i is maximal at ei < < 62. Expansion ([lOD is valid for any first order phase 
transition, as long as the correlation length ^ < L. At the critical temperature ^ remains 
finite, so that for sufficiently large L the appearance of a free energy gap AF indicates a 
first order transition; for small L < ^ the free energy is dominated by the bulk term. As the 
system approaches a tricritical point, ^ grows and the double minima structure can only be 
seen for large L. At a tricritical point and beyond it, the phase transition is continuous and 
hence there is no double minimum structure for any L. 



B. Results 

To obtain good statistics, we took histogram data every lO*'' MCS for 6 x 10^ MCS, 
averaging over 16 realizations of the disorder. 

Fig. 1^ shows two histograms of the internal energy density for p = 1 after normalization 
for two different temperatures (full and dashed line). The dots were produced by transform- 
ing the T = 1.77 data set (dashed line) to the lower temperature T = 1.74. The distribution 
depends sensitively on the chosen temperature; the Monte Carlo data at T = 1.74 and the 
transformed data from T = 1.77 show good agreement, indicating, that the equilibrium 
distribution has been well estimated. The simulation temperature in Fig. ^ was chosen very 
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close to the critical temperature. Therefore, as noted above, the distribution clearly has two 
maxima at ei and 62, which are equivalently minima of F{e). 

If we define the effective critical temperature Tc{L) by the equality of the two maxima 
-Ps(ei) and -P/3(e2) and extrapolate Tc(L 00), we obtain Tc{p = 1.0) = 1.716(5). This is 
slightly lower than the values from Ref. ^ (T^ = 1.736(1)) and Ref. ^ (Tc = 1.746(5)). This 
discrepancy may be due to using different definitions of "finite-size Tc" in the respective 
references. 

While for p = 1 the energy gap appears already at L = 10, for stronger dilution it 
can only be seen at larger system sizes, indicating the growing correlation length as one 
approaches the tricritical point. If a gap appeared at one system size, it continued to be 
present and in fact grew for larger sizes. The scaling behaviour according to eq. ([T0|), i. e., 
AF{L) oc L^, is fulfilled within the error margins. In the diluted case we found considerable 
fluctuations of the energy gap size depending on the realization. At p ~ 0.9 for our largest 
systems {L = 22) certain realizations were found, that showed none of the typical two peak 
structure, while others still showed a small gap. 

In order to determine the tricritical dilution ptrij we evaluated the percentage of real- 
izations with gap (= Y) with decreasing concentration of spins. For our largest systems 
(L = 22), at p = 0.91 Y = 90% of the realizations still showed a double peak, at p = 0.9 
only 40%, at p = 0.89 10% and for p = 0.88 no double peaks were found at all, i. e., 
all systems showed Gaussian peaks at T^. From extrapolating (Tc{L) , ptri{L)) determined 
by the condition Y = 0.5 for different system sizes (see dots in Fig. |10D, we estimated 
Ptri = 0.85 ± 0.03 in accordance with Ref. |19|. To determine the triciritcal concentration 
more exactly would require a precise scaling analysis based on more extensive data, which 
was beyond the scope of this work. 
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V. CONTINUOUS ANTIFERROMAGNETIC TRANSITION (P = 0.8) 



In this section we concentrate on simulations performed for p = 0.8. Our aim is to check, 
whether at the system orders into an antiferromagnetic state and, since we are below 
the tricritical concentration we measure critical exponents via finite size scaling close to the 
continuous phase transition, that we encounter. 



A. Quantities analysed 

For the antiferromagnetic phase, the staggered magnetization "vector" (see eq. H) is 
the appropriate order parameter. Thus, we calculate the second moment 

= [{m'^ ■ m^)]. (11) 

For T > Tc this is proportional to the staggered susceptibility x'' = Nm'^, which we analysed 
by using the finite size scaling form 

X\L,T)=L'-'^X\L'/''{T-T,)) (12) 

to extract the critical exponents t] and u and the critical temperature Tc. 

Since we are mainly interested in the magnetic order of the different dilution regimes, 
we also calculated correlation functions. To save computer time, we Fourier transform the 
lattices to fc-space and calculate the Fourier transformed correlation function, i. e., 

Gik) = [(Isfcp)] = l^exp(-in,-fc) (13) 

We applied the Fast Fourier algorithm to the (most relevant) L = 16 systems only and 
calculated the correlation function along the three (100) directions. For an antiferromagnet, 
the three k = n modes of G should yield the static staggered susceptibility, which we used 
as a consistency check. The correlation length ^ can be extracted from the knowledge of the 
scaling form of G, i. e., 

G{k) = ^^G{kO. (14) 
13 



The scaling factor G has the following asymptotics: For > and T T+, G{k^) — * 1 and 
and for T > Tc and A; — (in the afm case k = tc — k' with the limit k' — > n) 

G{k ^ 0) ^{kO'-"^ = (15) 

In eq. (|T5|) we have used the scaling form of the correlation length ^ ~ (T — Tc)^'^ and the 
scaling law 7 = z/(2 — r^). With this knowledge of the asymptotic form of G we chose the 
Ansatz 

that we used to fit our data of the correlation function in order to obtain an estimate for ^ 
and rj. The constant A has been introduced as the amplitude of the correlation function. 
Finally we shall analyse the heat capacity 



{e') - (e)' . (17) 



This quantity indicates a continuous phase transition to an ordered state by a weak diver- 
gence at the critical temperature. We use the finite size scaling form 

C(L, T) = L"/'^C'(Li/'^(T - T,)) (18) 

to extract critical exponents a and u as well as Tc. 

B. Results 

Our data of the staggered magnetization show an increasing (eq. |1T]) as T approaches 
the critical temperature, becoming more pronounced for larger lattice sizes. This contri- 
bution to the afm order parameter can as well be seen in the divergent behaviour of the 
staggered susceptibility (Fig. The scaling analysisiiof the staggered susceptibility yields 

14 



Tc{p = 0.8) = 1.07 (0.05), u = 0.51 (0.1) and rj = 0.05 (0.1). The data scale well over a wide 
temperature range with the exception of the L = 4 data and the errors in the exponents are 
within acceptable limits. 

The antiferromagnetic correlation length ^afm is found to grow continuously with de- 
creasing temperature and reaches half the lattice size at T ^ 1.2, before the critical temper- 
ature Tc = 1.07. Since finite size effects become appreciable at distances close to half the 
lattice size, we only included data above T = 1.2 for a fit of the scaling form ^ ~ (T — Tc)^'^. 
We extracted u = 0.55 by regression, which is in good agreement with our previous result. 



From fitting the correlation function to eq. (piq) , we were also able to extract a second 
estimate of rj; here, we found rj = —0.04 with a slow drift to ?7 — —0.02 for T 1.2. 

We do not perform a complete finite size scaling analysis of the correlation function in 
this work, because it would require a substantial amount of additional data for larger lat- 
tice sizes. Therefore, our estimates for the exponents obtained from G{k) are less reliable 
than those obtained from finite size scaling. Nevertheless, this analysis serves as an addi- 
tional consistency check and the critical exponents lie well within the error margins of those 
exponents obtained via finite size scaling of the susceptibility. 

The heat capacity shows a weak divergence at the critical temperature at p = 0.8. Fig. |^ 



displays the best result of the scaling procedure according to eq. (|T8|) with Tc = 1.08 (0.05) 



z/ = 0.57 (0.1) and a = 0.38 (0.15). These exponent values satisfy the hyperscaling relation 

du = 2-a (19) 

with u = 0.54 in 3 dimensions. 

In summary, the dilution regime p = 0.8 exhibits a continuous phase transition to an 
antiferromagnetically ordered state. We obtained critical exponents using finite size scaling; 
however, since p = 0.8 is very close to the tricritical point, it is quite possible that these 
are only effective exponents from the crossover between the mean field exponents of the 
tricritical point (?7 = 0, u = 1/2 and a = 1/2) to whatever universality class is appropriate 
for the continuous ordering transition. With the data at hand this question has to remain 
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open. 



VI. SPIN GLASS ORDER 

Upon further dilution we encounter a dramatic change in the magnetic order of the 
system at the phase transition. We investigate the question, whether below p = 0.8 the 
system really orders into a spin glass or if antiferromagnetic ordering can still be found. 
First, we introduce the spin glass order parameter and the spin glass susceptibility. Then 
we discuss our data at p = 0.7, and proceed with results from simulations of stronger diluted 
systems. 



A. Theory and quantities measured 

Besides the quantities that have been introduced in the previous section, we also mea- 
sured the standard indicator of a spin glass transition, the spin glass susceptibility xsg (see 
eq. @). For an infinite system a spin glass transition is signalled by a divergence of xsG 
as (T — Tc)~'^, with 7 = (2 — ri)^. For our scaling analysis, xsg is computed as the second 
moment of the overlap, defined by eq. @ and eq. @. We analysed our estimate of the spin 
glass susceptibility by using it's finite size scaling form 

Xsg{L, T) = L'~^ XsciL'/'iT - T,)). (20) 

Another important quantity, that is well-known in the analysis of spin glass simulation 
data0, is the Binder cumulant@ of the spin glass order parameter 

,.if3-MlV (21) 

It has the pleasant finite size scaling form 

g{L,T) = ~g{L'/''{T-T,)) (22) 

with no power of L multiplying g, which makes it very valuable for precise scaling analysis. 
The Binder cumulant (|2l| ) is defined so that < < 1, and above Tc, g{L,T) for 
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L oo. In particular, the intersection of all g{L,T) curves at some point provides an 
accurate estimate of T^. 

To investigate a change of magnetic order further, we analysed the correlation function 
again, this time also the spin glass correlation function 

Gscik) = [{\qk\')] = ^Eexp(-in, -fe) [(g.g,)]. (23) 

Here, we defined Qi = s[^\t + to)s'^\t + to)- The same fitting procedure of our data was 
employed as in the previous section in order to compute the respective correlation length, u 
and 77. 



B. Results for intermediate dilution (p = 0.7) 

The (antiferromagnetic) staggered susceptibility as shown in Fig. ^, is drastically reduced 
in comparison to Fig. H, and shows only a small tendency to increase as T decreases. 
Furthermore, scaling according to eq. (|12D could not be achieved for reasonable parameters. 

On the other hand our data reveal a divergence of xsg at p = 0.7 of the same order 
of magnitude as for p = 0.8 in Fig. ^. This divergence becomes particularly strong for 
larger lattice sizes in the vicinity of T = 0.85. If xsg is the critical quantity for this system 



then it should also scale according to eq. ED 



A finite size scaling analysis of the spin glass susceptibility is given in Fig. |^ All data 
for xsG with the exception of L = 4 scale well with Tc{p = 0.7) = 0.85 (0.05), z/ = 1.0 (0.2) 
and 7] = 0.1 (0.2). The fact that the spin glass susceptibility satisfies the above scaling form 
(PPP strongly suggests that the magnetic order of this model has changed between p = 0.8 
and p = 0.7 from antiferromagnet to spin glass. 

In Fig. 1^ we present our data of g at p = 0.7. The data show the typical behaviour 
as it has been observed in short range spin glasses, i. e., the data merge at approximately 
T = 0.85, indicating the phase transition. There is even a slight tendency of fanning out of 
the data below this intersection point, which is a strong evidence for the occurrence of a phase 
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transition. Such a fanning out is usually observed in uniform systems. In our system this may 
be due to the proximity to the tricritical point or to residual short range antiferromagnetic 
order (see below). The best fit was achieved with Tc{p = 0.7) = 0.83 (0.05) and z/ = 
1.05 (0.2), which agrees well with our estimates from the scaled spin glass susceptibility. The 
stronger scattering of this quantity, especially of the [{q*)] data, may be due to insufficient 
disorder averaging (40 systems for L = 16) and could probably be decreased with additional 
computational power. 

In Fig. D the temperature dependence of ^sg and ^afm are presented. While the spin 
glass correlation length is very small for large T but increases drastically as T ^ Tc ~ 0.83, 
the antiferromagnetic correlation length starts at a higher level but increases slower than 
^sg! it ceases to increase at about T = 0.9 and saturates for lower temperatures. 

All these results of our simulations at p = 0.7 are consistent with the interpretation 
that in this dilution regime we witness indeed a spin glass transition, but also encounter 
antiferromagnetic order of long but finite range being embedded into long range spin glass 
order. Apparently, a change of magnetic order has taken place within the interval p^, G 
(0.7,0.8), where denotes the critical dilution, where this change happens. 

C. Results for strong dilution 

Additional simulations were performed for concentrations p = {0.6,0.5,0.4,0.3}. With 
stronger dilution we find that the spin glass phase already encountered for p = 0.7 persists. 
Both, spin glass susceptibility and Binder cumulant of the order parameter scale well in this 
regime with slight dilution dependent critical exponents. Also, in this regime the data of 
the Binder cumulant stay together below Tc for all sizes, as has been observed in other short 
range Ising spin glass simulations. Unfortunately, it is increasingly difficult in this regime to 
equilibrate the systems. As the concentration is lowered, the critical temperature decreases 
rapidly (as expected from our argument earlier, that Tc{p'J = with p'^ > Pc = 0.195), but 
the characteristic microscopic energy barriers remain of order unity. Metastability becomes 
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an increasingly important problem and is prohibitive in large systems. In our simulations, 
for p < 0.6 we were unable to equilibrate the L = 16 systems sufficiently close to within 
reasonable computer time. 

From analysing the correlation function in the strong dilution regime we find two results: 
First, the spin glass correlation length shows a similar divergence close to the critical temper- 
ature for p = {0.6, 0.5, 0.4} as we have seen in the previous subsection, confirming again the 
development of long range spin glass order. Second, we still find short range antiferromag- 
netic order, which decreases with lower concentration: for p = 0.6 ^afm rises slowly when T 
is lowered and at Tc we have .^afm ~ 4; for p = 0.4 the correlation length remains constant 
at ^AFM ~ 2 for all temperatures which we can simulate (T G [0.56, 1.0], Tc{p = 0.4) ^ 0.47). 

Furthermore, both critical exponents u and rj apparently decrease with increasing dilution 
(see the fitted values in table I). The generic theoretical expectation is that the exponents 
should be universal all along the spin glass transition line, but the issue of universality is not 
completely settled in diluted systems^. In our case, we note that the drift in the exponents 
is pronounced for 0.5 < p < 0.7 but minimal for 0.5 > p > 0.3. Thus we suggest that 
this dependence is an artifact of the antiferromagnetic correlation length, which is large and 
p-dependent for 0.5 < p < 0.7, but is quite small for 0.3 < p < 0.5. However, we have 
insufficient information to check this proposition, and shall come back to this issue in the 
conclusion. 

We conclude, that for the whole range p G [0.3, 0.7] our simulations testify the existence 
of a spin glass phase transition for the short range antiferromagnetic Ising model on a fee 
lattice. Simultaneously, the model also exhibits residual antiferromagnetic order, that is 
relatively long ranged at p = 0.7 and saturates close to the critical temperature; it decreases 
with further dilution, showing no temperature dependence for p = 0.4. 
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VII. CONCLUSION 



In this paper we have presented a Monte Carlo simulation of the diluted short range 
antiferromagnetic Ising model on a fee lattice. The main purpose was to investigate the 
thermodynamic equilibrium properties of this model. 

In the undiluted case {p = 1) we have found a first order phase transition to an antifer- 
romagnetically ordered state, consistent with earlier simulations. Upon slight dilution, the 
antiferromagnetic order persists, the first order transition becoming weaker and changing to 
a continuous transition at the tricritical concentration ptri ~ 0.85. At p = 0.8 we still find 
antiferromagnetic order from analysis of the correlation function and scaling of the stag- 
gered susceptibility. Together with scaling of the heat capacity, we find mutually consistent 
critical exponents. The question of universality remains unclear; in view of the proximity 
to ptri ~ 0.85 it seems likely that we observe tricritical exponents at this point. 

Below p = 0.8 the divergence of the spin glass correlation length as well as scaling of 
the spin glass susceptibility and the Binder cumulant of the order parameter signal spin 
glass order. This means, that there must be a multicritical point at some concentration 
p* G (0.7,0.8), where the change of magnetic order takes place. Simultaneously, the quasi 
temperature independent staggered magnetization and the saturation of the antiferromag- 
netic correlation length suggest the breakdown of antiferromagnetic long range order. The 
fact, that the antiferromagnetic correlation length saturates while spin glass correlations 
still grow contradicts the view of a dynamically inhibited transition to an antiferromagnet. 
Rather, our data suggest that the coexistence of antiferromagnetic short range order to- 
gether with long range established spin glass order seems to be a special phenomena of this 
diluted model. 

Let us now turn to the question of universality in the spin glass phase. In this phase (see 
Sec. |V^ ) our scaling fits yielded critical exponents slightly dependent on dilution (compare 
table D; we speculated, that this could be an artifact of the changing antiferromagnetic 
correlation length. Even in unfrustrated models, there is still controversy over the univer- 
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sality of exponents under dilutionEl. Although our simulations are less extensive than the 
MC-simulations of the short range Ising spin glass by Bhatt and Young and by Ogielski, 
it is interesting to compare our exponents with their values for the ±J-model in ci = 3. 
For this model, they find0@ v ^ 1.2 and r] ^ —0.25. These values are just at or slightly 
out of the error margins of the present simulations, so that it is not entirely clear, whether 
the two models lie in the same universality class. Note also, that our value of 7 is smaller 
than for the "classic" spin glass simulations, roughly 7 ^ 1.8 for p G [0.3,0.7]. Additional 
simulations of larger systems would be desirable to confirm our result. 

It should be noted that even the best current simulations on hypercubic lattices have not 
put to rest the basic question whether the lower critical dimension for Ising spin glasses is 
below 3, or equal to 30@0iiilii. At the more modest scale of our simulation, it cannot be 
decisively answered whether the diluted fee Ising antiferromagnet has a genuine spin glass 
transition, and if so whether it has precisely the same exponents as the standard example 
of ±J spin glasses on hypercubic lattices using array processors. However, our data is 
consistent with both of these propositions. 

Although in this work we have performed an Ising model simulation we find striking 
similarities with experimental observations in DMS, which have Heisenberg like magnetic 
moments^. Not only are the respective phase diagrams (compare Fig. |lO| and Ref. H) 
qualitatively similar and the various critical concentrations are numerically close, but also 
the concentration dependence of ^afm in the spin glass phase agrees well with experiment. 
To be more specific, a comparison of our simulations with intensive experimental studies by 
Giebultowicz et al.i on Zni_pMnpTe shows the following agreements: (i) the afm transition 
is first order for p G [0.85, 1.0], (ii) the transition is continuous (to an afm long-range state) 
for p G [0.75,0.85], (iii) the afm order is of large but finite range for p < p^, ^ 0.75 in 
the spin glass phase and (iv) the afm correlation length decreases with increasing dilution 
p < p^. These results together with the measured divergence of the nonlinear susceptibility 
in Cdi_pMnpTe (see Ref. ^) below p^ ~ 0.75 do support the view of a spin glass phase in 
experimental Heisenberg systems. 
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TABLES 

TABLE I. Critical quantities for different concentration p from scaling of the respective sus- 
ceptibility, Binder cumulant and analysis of the respective correlation function. Not listed are the 
scaling results of the heat capacity for p = 0.8. They are: Tc = 1.07(0.05), a = 0.38(0.1) and 
V = 0.57(0.1). With susceptibility and correlation function we denote the valid quantity for the 
respective dilution regime. 

Critical Temperatures and Exponents 



Order Concentration Susceptibility Binder cumulant Correlation function 







Tc= 1.07 


(0.05) 








AFM 


0.8 


u=0.51 

//=0.05 


(0.10) 

(0.15) 






1] = -0.01 (0.05) 






rc=0.83 


(0.05) 


Tc=0.83 


(0.05) 




SO 


0.7 


u=im 

^ n 1 n 

f] — U.iU 


(0.20) 
[yj.zu ) 


u=1.05 


(0.20) 


u=0.96 

Tj — U.IO [U.ZU) 






Tc=0.75 


(0.05) 


Tc=0.76 


(0.05) 




SO 


0.6 


u=0.80 
r/=0.00 


(0.20) 
(0.20) 


z/=0.90 


(0.20) 


i/=0.98 
?7=0.07 (0.15) 






Tc=0.55 


(0.05) 


Tc=0.53 


(0.05) 




SO 


0.5 


u =0.73 
r] = -0.30 


(0.20) 
(0.25) 


u =0.73 


(0.20) 








rc=0.47 


(0.10) 


Tc=0.44 


(0.10) 




SO 


0.4 


u=0.80 
r] = -0.35 


(0.20) 
(0.20) 


u=0.80 


(0.20) 








Tc =0.27 


(0.10) 


Tc =0.27 


(0.10) 




SO 


0.3 


u =0.70 
r) = -0.40 


(0.20) 
(0.20) 


u =0.70 
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FIGURES 



FIG. 1. Internal energy density e versus temperature for p = 1. 

FIG. 2. Internal energy density versus temperature for p = 0.9. 

FIG. 3. The lines show two energy density distributions -P(e) versus energy density e for 
MC-simulations at two different temperatures, p = 1.0 and L = 16. The left peak corresponds to an 
antiferromagnetically ordered phase (low energy), the right peak to a paramagnetic phase, divided 
by an energy gap (latent heat). The points are calculated according to the Ferrenberg-Swendsen 
method from the high temperature distribution and agree well with the full line from low temper- 
ature MC-simulation. 

FIG. 4. Staggered susceptibility versus temperature, p = 0.8. 

FIG. 5. Scaled heat capacity CL~°'/^ versus L^/^{T — Tc), p = 0.8. Best scaling was achieved 
with Tc = 1.08 (0.05), u = 0.57 (0.1) and a = 0.38 (0.1). 

FIG. 6. versus temperature, p = 0.7. 

FIG. 7. xsgL^""^ versus lV'^(T -Tc),p = 0.7. Best scaling was achieved with Tc = 0.83(0.05), 
1^ = 1.0 (0.2) and r? = 0.1 (0.2). 

FIG. 8. Binder cumulant of the spin glass order parameter versus temperature, p = 0.7. 

FIG. 9. Spin glass correlation length ^SG (*) and antiferromagnetic correlation length ^afm 
(•) versus temperature, p = 0.7. 
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FIG. 10. Phase diagram of the short range fee Ising antiferromagnet with dilution. We plot 
Tc versus spin eoneentration p. We observed a I''* order phase transition to an afm-state for 
1.0 > p > ptri ~ 0.85 (long dashed line), becoming continuous for ptri > p > p* ^ 0.75 (dotted 
line). The dots above the transition line denote points {Tc{L) , ptri{L)) with the condition Y = 50% 
(see Sec. IV). Upon further dilution {p < p*) we encountered a spin glass transition (full line, being 
extended down to the percolation threshold pc = 0.195). The lines serves only to guide the eye. 
The gothic arch marks a region, where the order of the system is unknown (afm, spin glass or 
coexistence of both) . 
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